Theory of magnetic field-induced metaelectric critical end point in BiMn 2 Os 
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A recent experiment on the multiferroic BiM^Os compound under a strong applied magnetic field 
revealed a rich phase diagram driven by the coupling of magnetic and charge (dipolar) degrees of 
freedom. Based on the exchange-striction mechanism, we propose here a theoretical model with the 
intent to capture the interplay of the spin and dipolar moments in the presence of a magnetic field 
in BiMn2 05. Experimentally observed behavior of the dielectric constants, magnetic susceptibility, 
and the polarization is, for the most part, reproduced by our model. The critical behavior observed 
near the polarization reversal (P = 0) point in the phase diagram is interpreted as arising from the 
proximity to the critical end point. 
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I. INTRODUCTION 

Beginning with the pioneering work of Hur et al. , a 
series of experiments has uncovered remarkable cross- 
correlations of the magnetic and electric dipole (i.e. po- 
larization) behavior in a class of compounds RMn 2 05 
(R=Tb,Ho,Dy) 2 ~ 4 . The coupled behavior of the mag- 
netic and polarization degrees of freedom is due in large 
part to a significant exchange-striction in these mate- 
rials, and to the presence of geometric frustration in 
the magnetic exchange network. The idea of exchange- 
striction as the driving force of multiferroic behavior in 
the RMn 2 05 compound was proposed in Ref. 3. 

A recent high magnetic (H) field study on one mem- 
ber of the RMn 2 5 family, BiMn 2 5 (BMO), revealed a 
high-field phase with critical behaviors of the polarization 
and the magnetization at the point where P (bulk polar- 
ization) is tuned through zero 5 . In the low-temperature 
ferroelectric phase of BMO, application of the magnetic 
field H along the crystallographic a-axis in excess of 20 
Tesla resulted in a sharp increase in the &-axis dielectric 
constant, as well as in the slope of a-axis uniform magne- 
tization dM/dH, as the field strength swept through the 
critical value H c . The temperature(T)-dependent trace 
H c (T) agreed well with the position of P — separating 
the low- field P > from the high-field P < region 6 , as- 
suming that the H = state had the P > polarization 
to begin with. Down to the lowest temperature measured 
at 0.66K, the P > to P < crossover appeared to be 
smooth with no sign of a first-order discontinuity. Fur- 
thermore, the behavior of P at 0.66K near H = H c was 
shown to agree well with the power-law \P\ ~ \H— H c \ x / 3 , 
while that of the 6-axis dielectric constant was repro- 
duced with e b (H)-e b (H = 0) - \H-H C \~ 2 ^. A 
Ginzburg-Landau scheme was employed to explain the 
observed power-law behavior 5 . 

As is obvious from the symmetry consideration, a 
second-order phase transition at P = is ruled out be- 
cause both sides of P = are already symmetry-broken 
states. Only a first-order discontinuity or a crossover is 



left as a possibility. It was then conjectured 5 that a criti- 
cal end point with an extremely low critical temperature 
T* must exist in this material. The observed critical be- 
havior in both P and £& at low temperature then follows 
naturally from the proximity to the putative critical end 
point, it was claimed 5 . 

Given the novelty of the claim and excitement over the 
possible metaelectric phenomena in a multiferroic com- 
pound, it is desirable to develop a microscopic model that 
can capture the essential aspect of the observed dielectric 
and magnetic behavior of BMO under a high magnetic 
field. While the model we propose is based on the exist- 
ing exchange-striction ideas of Refs. 3 and 5, this is the 
first attempt to examine the exchange-striction physics in 
RMn 2 05 at a microscopic level. In Sec. II, the complex 
structure of magnetic Mn networks for BMO is reduced 
to a simple, manageable spin model coupled to lattice 
displacements. The model naturally embodies the ideas 
of spin-lattice coupling already proposed for other com- 
pounds such as YMn 2 05 3 . The relation of the frustration 
in the magnetic exchange network to the local displace- 
ment of Mn ions is made transparent. Then in Sec. Ill a 
thorough classical Monte Carlo simulation of our model 
is carried out, both justifying the continuous spin flop 
model introduced in Ref. 5 and revealing the power- 
law behaviors of susceptibilities as in the experiment. 
The observed exponents agree fairly well with the ex- 
perimentally measured values even though no quantum- 
mechanical consideration is given in the present model. 
The phase diagram for our model is indeed consistent 
with the presence of a critical end point. We close with 
a summary and outlook in Sec. IV. 

II. THE MODEL 

The pronounced feature of the magnetic structure of 
BMO is the geometrically frustrated nature of the mag- 
netic interaction pathways. The Mn atoms in BMO oc- 
cur in two varieties: Mn 3+ (whose spin is S = 2 and is 
surrounded by an oxygen tetrahedron) and Mn 4+ (spin 



FIG. 1: (color online) Network of Mn atoms in BiM^Os. 
Filled and empty atoms are Mn 3+ (5" = 2) and Mn 4+ (S = 
3/2), respectively. A unit cell containing eight Mn atoms is 
shown as a cube with its axes labeled a, b, and c. Four unit 
cells are shown in the figure. Bars connecting the atoms have 
non-zero exchange energies. Exchange interaction between 
the two Mn 4+ atoms will be ignored. That makes the unit 
cell with six independent spins. 




FIG. 2: (color online) Projection of the Mn network onto the 
ab plane with six atoms per unit cell. Thick and thin full lines 
represent J5 and J4 bonds, while the green dotted lines are J3 
bonds. A sample spin configuration with R(ight) and L(eft) 
pointing spins are displayed. The J3 bonds alternate between 
being fully satisfied and fully frustrated. 

S = 3/2, surrounded by an oxygen octahedron). The 
large spins of both Mn atoms allows us to treat them as 
classical to the first approximation. 

The real-space locations of Mn atoms and their ex- 
change network is presented in Fig. 1. There are eight 
Mn atoms in a unit cell with four Mn 3+ and four Mn 4+ 
ions each. Three antiferromagnetic exchange interactions 
have been identified in the literature as dominating the 
magnetic structure 2 . The two adjacent Mn 3+ ions (filled 
circles in Fig. 1) form the strongest exchange bond with 
J5. The exchange interaction involving one Mn 3+ and 
one Mn 4+ lying adjacent to it along the a-axis is the 
next strongest with J4. Magnetic exchange of Mn 3+ with 
Mn 4+ lying along the &-axis is given by J3, which is the 
weakest. All three J's are antiferromagnetic. As seen 
in Fig. 1, a given Mn 3+ spin is exchange-coupled to an- 
other Mn 3+ spin on one side (J5), and a pair of Mn 4+ 
spins on the other (J 4 ). The two Mn 4+ spins interact 
only weakly, and we will ignore this weak exchange of 







FIG. 3: (color online) A schematic representation of the Mn 
network. Thick and thin horizontal links are J5 and J4 bonds. 
The inter-chain bond J3 is shown as dotted lines. Two types of 
alternating chains are labeled as A and B. A unit cell contains 
six spins labeled 1 through 6. The magnetic unit cell is twice 
as large (shaded region) . Two kinds of Mn 3+ pairs, formed by 
2 — 3 and 5 — 6 atoms, exist in a unit cell. In the experiment 
of Ref. 5, a magnetic field is applied along the a-axis as 
shown and polarization develops along the 6-axis. The spin 
orientations are antiferromagnetic within a chain, and point 
in the direction dictated by the local anisotropy, which are 
different for the two chains. 



Mn 4+ spins for the sake of simplicity. As a result, the 
two Mn 4+ spin behave identically and there are only six 
independent spin degrees of freedom in a unit cell. The 
approximation to keep only J3, J4, and J5 also makes the 
system two-dimensional. 

The six independent spins in a unit cell are coupled to 
one another in the manner depicted in Fig. 2, where a 
zigzag chain consisting of alternating J5 — J4 — J4 — J5 — 
J a — J a — ■ ■ ■ bonds is shown running along the a axis. An 
antiferromagnetic spin configuration is realized for each 
chain. A weak antiferromagnetic coupling J3 exists be- 
tween the chains for a selection of Mn sites connected by 
dashed lines in Fig. 2. The situation is further simplified 
in the schematic plot of Fig. 3. Here the geometrically 
frustrated nature of the Mn exchange is apparent in the 
form of a closed loop consisting of five Mn spins. Because 
of this unique connectivity, the inter-chain interaction 
cannot be fully satisfied for all J3 bonds. For a particu- 
lar realization of antiferromagnetic order on the chains, 
the inter-chain antiferromagnetic interaction is alterna- 
tively fully satisfied and fully frustrated as one can see in 
the sample spin configuration of Fig. 2. Translating the 
spin configuration by one atom for a given chain merely 
shifts the locations of the frustrated bonds by one lattice 
atom, but fails to relieve the frustration itself. And as 
a consequence of the frustration, the ground state would 
possess 2 N degeneracy, N being the number of chains. 

In BMO as in other RMn 2 05 compounds, the frustra- 
tion is relieved through the spin-lattice interaction. For 
a given Mn 3+ pair (a pair of adjacent Mn 3+ ions), one 
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Mn 3+ spin is favorably exchange-coupled (anti-parallel 
spins) with the Mn + spin connected to it, but the other 
Mn 3+ spin must be unfavorably coupled (parallel spins) 
with its neighboring Mn 4+ spin. Then the Mn 3+ pair 
as a whole moves in the direction that strengthens the 
favorable bond. The relative positions of the Mn 3+ ions 
within a pair is assumed to remain rigid during the dis- 
placement, while the center-of-mass of the pair is allowed 
to move. If all Mn 3+ pairs are displaced in the same di- 
rection, one has a net polarization and a ferroelectric 
state. There are two types of Mn 3+ pairs in a unit cell, 
namely 2 — 3 and 5 — 6 pairs in Fig. 3. Although their 
movements are not strictly along the b axis in the real 
compound, it is also known that the a component of the 
displacements cancels out between the two Mn 3+ pairs, 
leaving only the b component to manifest itself in net 
polarization 2 . In this regard, BMO behaves as a uniaxial 
ferroelectric. 

The unit cell contains six independent spin sites la- 
beled 1 through 6 in Fig. 3. The spin-spin interaction 
energies within the chain (Ei) and between the chains 
(E2) read, respectively, 

El = J5 y^(Sj2 -Sj3 + Sj5-Sj6) 

i 

i 

E2 = J3 y^(Sji ■Sa 6 + S i 2-S i 4 + S i 4-S^ T 3 + Si5-S^ ! i), 

(2.1) 

repeated over all unit cell index i. Adjacent cells along 
the a- and 6-axes are labeled i±x and i±y, respectively. 

The spin-lattice interaction ties the displacement of 
the Mn 3+ pairs, or the local dipole moment, with the 
Mn spin configurations. Each unit cell i contains two 
Mn 3+ pairs. The displacement of the 2 — 3 and 5 — 6 
pairs along the fe-axis, labeled as di and Ui, are subject 
to the force generated through exchange-striction. There 
is also a potential energy increase associated with the 
displacements that, up to fourth order, can be written as 

£<(«? + d f)/^X + (7/4) Ei(«i + 4), where X P^ys the 
role of bare dielectric susceptibility and 7 is the interac- 
tion strength. With the suitable re-definition of x, Ui, di, 
and 7, one can define the strength of the spin-lattice cou- 
pling to be one, and arrive at the spin-lattice interaction 
energy 

i i 

— di(Si3 ■ Si-yA — Si2 ■ Sj4) 

i 

— ^ Wj(Sji • S i6 — S,5 • Sj^i). (2.2) 

i 

The last two lines express the exchange-striction effects. 
Because of the rescaling, we can regard \ both as the 



bare dielectric susceptibility and the spin-lattice coupling 
strength. 

To the above energies one adds the single-ion 
anisotropy contribution 

3 6 
Ei = -I^^(S, a . IJ i) 2 -J^^(S !a .™ B ) 2 . (2.3) 

i a—1 i a— A 

The local anisotropy axes ha and n b are assumed differ- 
ent for the A and B chains. Finally, one adds the Zeeman 
energy 

6 

E 5 = -Hj2J2 Sia ^- ( 2 - 4 ) 

i a—1 

The total energy governing the behavior of spins and dis- 
placements in BMO reads 

E = Ei+E 2 + E 3 + E 4 + E 5 . (2.5) 

This is the proposed "minimal model" for the BMO. In 
the subsequent section we do a classical Monte Carlo sim- 
ulation of this energy form. 

The bulk polarization P is due to the net displacement 
of the Mn 3+ pairs, 

P-^T^ + d,)- (2-6) 

i 

If we can ignore the quartic interactions in Ui and di, 
the dependence of the local displacements Ui and di on 
the surrounding spin configuration can be worked out 
exactly, and gives the polarization 

P OC ^ (S,3 • S^4 — S i2 • Sii) 

i 

+ ^ (S,i • S l6 — S i5 • S»tt/,i) • (2.7) 

i 

Before closing this section it is important to emphasize 
that the present model is purely classical in its nature. 
A proper quantum analogue will be worked out in the 
future. 



III. MONTE CARLO CALCULATION 

An antiferromagnet with the magnetic field applied 
along the direction of the single-ion anisotropy undergoes 
a spin-flop process at the critical field H c — \fJl, where 
J and I are the exchange and local anisotropy energies, 
respectively. If the field direction is not aligned with 
the anisotropy direction, the spin-flop occurs instead in 
a continuous manner as the spins gradually rotate with 
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H . Such a continuous spin flop can occur in BMO be- 
cause the local anisotropy directions ua and hs are not 
strictly parallel to the a axis, the direction of the applied 
field, but are off by ±8° 5 . The unique feature of BMO 
that follows from the different anisotropy directions of 
the two types of chains (A and B in Fig. 3) is that the 
spins on the two chains can rotate in the opposite direc- 
tions with increasing H. If indeed one set of chains has 
its spins rotate counterclockwise and the other set clock- 
wise, the once anti-parallel pair of spins becomes parallel 
and the parallel spins anti-parallel at sufficiently large 
field strength, and due to a relation such as Eq. (2.7), 
the polarization direction will get reversed. 

The salient features of the high-field experiment on 
BMO 5 are summarized here to facilitate the comparison 
with the Monte Carlo results. 

• The bulk polarization P along 6-axis reverses its 
direction at a critical field H = H c applied along 
the a-axis. Near P = and at the lowest measured 
temperature T — 0.66K, the field dependence of P 
is consistent with \P\ ~ |i7 — i? c | 1 / 3 . 

• The 6-axis dielectric constant ej, shows a pro- 
nounced peak as H is tuned through H c . The be- 
havior at T = 0.66K is consistent with £b(H) — 
e b (H = 0)~ \H-H C \- 2 I\ 

• The a-axis magnetic susceptibility also shows a 
peak at H = H c . 

• The temperature dependence of £&(T) with the field 
value fixed at H ss H c follows a non-Curie- Weiss 
form, known as the Barrett's formula 7 . 
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FIG. 4: (color online) Schematic H — T phase diagram of 
the model Eq. (2.5) for (a) weak, (b) moderate, and (c) 
strong spin-lattice coupling \. The dashed and full lines sep- 
arating the P > from P < ferroelectric (FE) region are 
crossover and first-order transition lines, respectively, and the 
dark square in (b) is the critical end point. The scenario (b) 
is most consistent with known facts about BMO. 

The full lattice model of Eq. (2.5) was treated within 
the classical Monte Carlo scheme to see if the above- 
mentioned features of the experiments can be captured 
within our model. Aided by the experimental input, we 
consider the planar spins confined in the ab plane, and 
work with the two-dimensional lattice disregarding the 
coupling along the c-axis. A lattice of L x x L y unit cells, 
each unit cell consisting of six Mn sites, is considered. 



We choose the field directed along the a axis as in the 
experiment 5 , and let H vary from to +H max for each 
fixed temperature. The calculation was then repeated for 
many different temperatures. -ff max is chosen in such a 
way that P evaluated from Eq. (2.6) or Eq. (2.7) van- 
ishes before \H\ = H max is reached. Such a field- induced 
paraelectric transition was continuous, and occurred be- 
fore the full polarization of spins due to the strong Zee- 
man field could take place. 

A difficulty with the present simulation is the lack 
of information about the parameter values such as J3 
through J5 and spin-lattice coupling strength \- Initially, 
we worked with several different sets of parameters and 
later identified the ones which best reproduce the ex- 
perimental facts. In the course of the general search, 
we realized that three distinct behaviors (Fig. 4) are 
possible for the P > to P < crossover: (a) With 
a sufficiently weak x, t ne entire P = line becomes 
a crossover without a discontinuous jump in P at any 
temperature, (b) The intermediate range of x gives the 
P = curve that begins as a first-order critical line at 
low temperature but terminates at a finite temperature, 
T* , as a critical end point. The higher-temperature part 
of the curve becomes a crossover, (c) For a sufficiently 
strong spin-lattice coupling x, the entire P = line is a 
first-order transition that merges with the second-order 
paraelectric transition line at high temperature. It is 
the behavior near the critical end point in scenario (b) 
that is most relevant for BMO. The Monte Carlo re- 
sults discussed below are for the parameters that give 
rise to the scenario (b): J4/J3 = J5/J3 = 20, X/J3 = 3, 
I/J3 — 9, and 7 = 0. The anisotropy angles 9a and 9b, 
defined by ua • x = cos 9a , Ub • % — cos 9 b , are chosen 
as 9 a = —9b = 30°. The exaggerated anisotropy angle 
(experimental values are ±8°) is a consequence of search- 
ing for a parameter set that can produce the critical end 
point temperature T* at a sufficiently low temperature, 
well below the paraelectric transition. 

All calculations were performed on the lattice size 
L x = L y = 16 with the periodic boundary conditions in 
both directions. A standard Metropolis update scheme 
was used. Due to the complexity of the model, some care 
was needed in implementing the Monte Carlo program. 
First, a "typical" configuration at each temperature T for 
zero field was obtained by means of simulated annealing 
method. Then, beginning with the zero-field configura- 
tion thus obtained, we increase H to compute physical 
quantities as functions of T and H. At each temperature 
and field, at least 2 x 10 4 Monte Carlo steps per spin and 
displacement were made, and typically 4 x 10 3 steps were 
discarded for equilibration. Near the critical region, more 
steps of up to 5 x 10 5 were required for sufficient equi- 
libration and ensemble averages. Throughout the paper 
we denote energy, temperature, and field in units of J3. 

For a given instantaneous configuration, we compute 
the magnetization per spin along the field, 
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FIG. 5: (color online) (a) Polarization P and (b) uniform o- 
axis magnetization M as a function of magnetic field H at 
various temperatures T. The critical end point occurs be- 
tween T = 0.1 and T = 0.2. 



and the polarization per unit cell 

V = J—lj2( u i + d i)- (3-2) 

The average polarization P and magnetization M is then 
calculated by 

P = (V),M = (M), (3.3) 

where (...) indicates the ensemble average. We can also 
compute the dielectric (xp) and magnetic (xm) suscep- 
tibilities as 

XM = 6 ±±f((M>) - (Mf). (3.4) 

Varying the parameters within the scenario (b) of Fig. 4 
only gave rise to minor quantitative differences without 
altering the main results described below. The reduc- 
tion of \6a\ — \&b\, for instance, resulted in the overall 



increase of \P\ and enhanced T*. The introduction of 
nonzero 7 only reduces \P\. For these reasons we believe 
the results presented in the following represent the gen- 
eral features near the critical end point in scenario (b). 

In Fig. 5 the polarization P is plotted against H for 
various temperatures. The behavior at T — 0.1 showed a 
jump from P > to P < as in a first-order transition. 
The corresponding a-axis magnetization also undergoes a 
sudden increase at H — H c . For T > 0.2, both M and P 
evolve continuously with a sharp slope at H = H c . The 
critical field position H c itself depends smoothly on the 
temperature. We note that H c deduced as the location of 
P = in the P vs. H plot is numerically slightly different 
from the positions of the maximum susceptibilities. The 
same difference also shows up in the experiment 5 , but 
we do not have a good reason to believe that the small 
discrepancy has any physical importance. 
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FIG. 6: (color online) (a) Dielectric susceptibility \p an d (b) 
uniform magnetic susceptibility \m as functions of magnetic 
field H and temperature T. The peak occurs near where 
P = 0. The peak height rises upon lowering the temperature. 
The lowest-temperature peak at T = 0.1 (just below T*) is 
smaller than the peak at T = 0.2 (just above T*). The second 
set of peaks at higher magnetic fields are due to the ferro- to 
para-electric transition. 

The susceptibilities xp an d Xm from Eq. (3.4) are 
shown in Fig. 6. Clear peaks in both quantities were 
found as H crosses H c , and the heights of both peaks 
increased upon approaching T* from above. Both are 
expected to diverge at the critical end point (H*,T*). 
The peaks grew smaller at T = 0.1, which lies below 
T* . In the experiment both susceptibilities reached max- 
imum peak heights at ~ 5K and decreased below it. On 
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the other hand, no sign of a first-order transition was 
found for temperatures below 5K, and no sign of diver- 
gent susceptibilities at or near 5K. Hence it is incorrect 
to conclude that ~ 5K corresponds to T* in the exper- 
iment. Rather, the genuine first-order transition should 
take place, if at all, below the currently available temper- 
ature of 0.66K. It may be that the decrease of the sus- 
ceptibility that begins with 5K is a quantum effect such 
as the presence of a localized phonon of finite energy. 

The polarization P and dielectric susceptibility \p a * 
T = 0.2 (just above T*) and in the vicinity of H = H c 
are further analyzed in Fig. 7. Displayed on a log-log 
plot, the data are consistent with the power-law expo- 
nents a' = 1/3 and 7' = 2/3, the same exponents used 
to fit the experimentally observed behavior of P and Sb 
at T = 0.66K. A Ginzburg-Landau argument predicting 
the same exponents can be found in Ref. 5. 
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FIG. 8: (color online) Dielectric susceptibility \p (in arbi- 
trary unit) as a function of temperature T for various fields. 
The lines are best fits to the Barrett formula with the two 
temperature scales To and Ti as fitting parameters. 



fit several curves of \p versus T to the same formula in 
the vicinity of H C (T = 0)/J 3 w 27.0, the critical field 
value at zero temperature. For H < H c (0) (lower panel), 
it is apparent that the Barrett formula does not describe 
the curves very well. For H > H c (0) (upper panel), the 
curves seem to fit reasonably well to the formula, only if 
we allow for negative values of Tq although Tq should play 
the role of the critical temperature of the ferroelectric 
transition and remain positive. In contrast, the fit to the 
experimental data were made with positive To in Ref. 5. 
Overall, we do not find a good agreement of our Monte 
Carlo data for xp to the Barrett formula. A Curie- Weiss 
fit to the high-tcmperature side of the data also resulted 
in negative To ■ To achieve improved agreements between 
theory and experiment in this regard, we believe it is 
essential to consider the quantum-mechanical nature of 
the phonon modes u$ and di. 



(b) 

FIG. 7: (color online) Polarization P and dielectric suscep- 
tibility xp as functions of magnetic field H at temperature 
T = 0.2. The dotted lines represent the power-law behaviors 
\P\ oc \H-H c \ a ' and X p oc \H-H c \-~*' with the critical field 
H c /J 3 = 26.92 and the exponents a' = 1/3, 7' = 2/3. The 
errors are at most twice as large as the symbol. 

The quantum nature of the displacive phonon mode is 
reflected in the modification of the Curie- Weiss behavior 
of the dielectric susceptibility to the one described by the 
Barrett's formula 7 : 



Xp(T) = 



M 



(Ti/2) coth(Ti/2T) - T 



(3.5) 



It was shown that the experimental data for £&(T) fit 
well to the above formula 5 . In Fig. 8, we attempted to 



IV. SUMMARY AND OUTLOOK 

In this paper, we proposed a minimal model of mag- 
netic field-induced metaelectric critical end point recently 
observed in B1M112O5. A classical energy involving the 
lattice and spin degrees of freedom and their coupling 
was written down in Eq. (2.5) and its properties an- 
alyzed with the Monte Carlo method. Our findings are 
summarized below. The readers will find it useful to com- 
pare the following set of results with the summary of the 
experimental facts given at the beginning of Sec. III. 

• The bulk polarization P along 6-axis did reverse its 
direction at a critical field H = H c applied along 
the a-axis. The spins on the A and B chains rotated 
continuously, and in the opposite directions, under 
the increasing field. Near P = 0, and at T slightly 
above T*, the field dependence of P was found to 
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be in reasonable agreement with the power-law be- 
havior, \P\ ~ liJ-i^l 1 / 3 . 

• The b- axis dielectric susceptibility \p shows a pro- 
nounced peak as H is tuned through H c . The be- 
havior at low temperature just above T* is consis- 
tent with xp ~ \H-H C \- 2 I 3 . 

• The a-axis magnetic susceptibility also shows a 
peak at H = H c which reaches a maximum value 
at T*. 

• The temperature dependence of Xp(T) at a fixed 
field H w iJ c is generally inconsistent with the 
Barrett's formula 7 . The experimentally observed 
Xp(T) agreed better with the Barrett's formula. 

In conclusion, the magnetic field dependence of the 
polarization, and magnetic and dielectric susceptibilities 
obtained from our model proved to capture most of 
the features of the experiment. The simultaneous rise 
in the dielectric and magnetic susceptibilities in the 
continuous spin flop regime emerges naturally from 
our model. Other features such as the temperature 
dependence of the dielectric susceptibility do not agree 
well with the experimental results. The height of the 
susceptibility peaks reaches a maximum at T* in our 
theory since that is where the expected divergence 



should take place, but, experimentally, the peak heights 
reach a maximum at ~ 5K without showing signs of 
a first-order transition below that temperature. These 
discrepancies calls for a refinement of the present model 
that should include, among other things, the quantum 
nature of the displacive phonon modes expressed as 
di and m in Eq. (2.5) and the quantum dynamics of 
the spins. To what extent the quantum correction will 
alter the low-temperature behavior of the classical result 
remains to be explored. It is encouraging, on the other 
hand, that a simple classical model such as we propose 
already captures many of the prominent features of the 
experiment. 
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